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Abstract 

In this note we study the zero temperature dynamics of the Sherrington Kirkp- 
patric model and we investigate the statistical properties of the configurations that 
are obtained in the large time limit. We find that the replica symmetry is broken 
(in a weak sense). We also present some general considerations on the synchronic 
approach to the off-equilibrium dynamics, which have motivated the present study. 



1 



1 Introduction 



In recent years there have been many progresses in our understanding of the non-equih- 
brium dynamics [l]-[6] of the infinite Sherrington Kirkpatrick model [7, 8] and of other 
glassy systems [9, 10, 11]. The aim of these studies is to compute the properties of the 
systems (e.g. the energy or the magnetization) as function of time if we know the initial 
configuration at time zero. 

These progresses have been done using a diachronic approach in which the evolution of 
the system is studied by writing closed equations for the correlation functions and response 
functions at different times. 

An alternative synchronic approach has been put forward [12, 13]; here we consider the 
probability distribution of the configurations of the system at a given time t (Pt(cr)) and 
we write it as 



In this approach there are two crucial steps: 

• (a) The determination of the effective Hamiltonian Ht- It is possible that in order to 
obtain qualitative and semiquantitative information it is not necessary to compute 
exactly Tit, and an approximate knowledge is sufficient. 

• (b) The computation of the statistical properties of the system at given effective 
Hamiltonian. 

This approach may be successfully only if the effective Hamiltonian Tit (or a reasonable 
approximation to it) is not too complicated. While the diachronic approach is rather 
systematic, a good amount of guesswork is needed in the synchronic approach in order to 
chose a reasonable form of the effective Hamiltonian. 

In order to get some intuition on the possible forms of the effective Hamiltonian we 
have studied in this note the statistical properties of the configurations that are obtained at 
large times in the SK model using a zero temperature dynamics. This problem has already 
been studied in the past [4, 6, 14, 15, 16, 17], and it has its own interest. We have obtained 
some new and unexpected results, i.e. we have found that the replica symmetry is broken 
(in a weak sense [18, 19]) and that the distribution of the local fields has uneexpected 
properties. 

In the second section we present some general considerations on the synchronic ap- 
proach. These considerations have been the motivation of the present study, but they may 
be skipped by the reader interested only in the results of the paper for the SK model. In 

the third section we define the dynamics which we study and we recall some known results 
on the zero temperature solutions of the TAP equations [20]. In section IV we present 
the numerical results of our investigations and finally in the last section we present some 
tentative conclusions. 



We consider a system whose variables satisfy some kind of stochastic or deterministic 
evolution equations. The probability distribution of the configurations at time zero is given: 



Pt{a) oc exp{-nt{cr)). 



(1) 



2 The Synchronous Approach 



Po oc exp(-7io)- 



(2) 
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The quantity TYq plays the role of boundary condition. 

In the simplest case the system is random at time zero and Hq = 0. Our aim is to 
compute Tit given TIq. If we are rather luckily we can compute it exactly. This happens 
for example in the Gaussian SK [13] or in the spherical SK model. Before discussing the 
general approach it may be useful to show in details the soluble example. 

2.1 The Gaussian Sherrington Kirkpatrick model 

In the Gaussian Sherrington Kirkpatrick [13] the Hamiltonian is 

H ^ -^Ji^kC^i ak, (3) 

i,k 

where the variables cxj are real, with an intrinsic probability distribution at infinite tem- 
perature equal to the product of uncorrelated Gaussian of unit variance. 

The variables Ji^k are usually Gaussian distributed and uncorrelated. The precise form 
of the distribution of variable J is not important here; we will only assume for simplicity 
of notations that the matrix J does not have zero eigenvalues. 

The dynamic is given by the Langevin equation: 

^^-a, + pF,{t)+v,it), (4) 



where the force is given by 



dH 

Ht)--^-Y.JrMi)k- (5) 



The variables r] are have an uncorrelated white noise distribution 

{Vi(t)Vk(t')) = 2T5{t - t%,k- (6) 

Wc suppose that at time zero the initial condition is simply a{0)i = 0. It is easy to 
show that the effective Hamiltonian at time t is given by 

nt{a) = J2aif{J)i,k<Jk, (7) 

i,k 

where the function /* is given by 

/'(^) - 13^^ (8) 

and f\J)i,k denote the i, k element of the matrix /*( J). 

The proof of this statement may be obtained by noticing that the evolution equations 
becomes much simpler in the basis where the matrix J is diagonal [13, 21, 22]. In this 
basis it is easy to sec that the components of the variables a along the eigenvectors of the 
matrix J are uncorrelated. Their variance can be computed and in this way one obtains 
equation (8). 
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2.2 The general approach 

In the general case the system has an Hamihonian Hj(a) and the evolution is described 
by a Langevin equation of the form in eq. (4) or by some analogous equation for discrete 
systems. 

We suppose that the effective Hamiltonian can be written as 

no^g{<j,J,X{t)), (9) 

where g is a preassigned function which depend on M variables Xa{t) {a — 1...M); M may 
also be infinite. In this framework the effective Hamiltonian depends on time only trough 
the variables X{t). 

If we suppose to know the function g (or a good approximation to it), the problem 
consists in finding the appropriate values of the functions Xa{t). 
There are two strategies which we can follow: 

• (a) We choose a set of observables Oa((7, J) for a = 1...M and we impose the validity 
of the equations: 

d{Oa)t _ d{Oa)t dX^ 

dt ~ ^ dX^ dt' ^ ^ 

where the r.h.s. is computed using the Langevin equatios. 
If we are new to equilibrium it is convenient to write 

Ht^PH + ^O^X^it). (11) 
For small values of Xj{t) we can linearize the equations and we get 

EiOaO,)/-^ = (12) 

where the expectation values are computed at equilibrium and ()c denotes the con- 
nected expectation value. 

• (b) We write down the Fokker Plank equation: 

f = £P(T). (13) 

where JC is the appropriate linear operator. A variational principle may be used to 
chose the variables Xn{t). For example we can impose that 

j2ij:^^-mt)r (m) 

{a} 7 "^^7 

takes the minimum value. 

In the past this point of view has been advocated in the study of turbulence and both 
strategies have been followed [23, 24]. 
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In their original work Coolcy and Sherrington [12] have analysed the dynamics of the 
SK model with Ising spins. They have taken M = 2. In the case of all spins equal to 1 at 
time zero, they have 

g{a, J, X{t)) = Xi{t) Ji,kcricrk + X2{t)Y^ai (15) 

i,k i 

Reasonable results have been obtained expecially at the short times. 

It is clear that the previous equation (with (M=2)) is only a first approximation. In 
the case of the Gaussian SK model it was proved [13] that it does not reproduce the exact 
result. Indeed in this last case we have seen that the correct expression (for (7i(0) = 0) is 
given by 

g{a,J,\{t))= Ut)Y.iJ"kk'^i'^k, (16) 

a=l,oo i,k 

where as usual {J")i,k is the matrix element of the matrix J to the power a. 

For the Ising case is possible that the previous formula is not adequate and that other 
terms may be present as 

(17) 

i k 

It is quite reasonable that near to the critical temperature too high powers of the 
variables a are not present so that it possible that a manageable approximation may be 
done using the techniques developed in [11]. 



3 Zero temperature dynamics and the TAP equations 

Here we consider the SK Ising model with the Hamiltonian in cq. (3), with the J random 
independent Gaussian variables with variance -j^, N being the total number of Ising spins. 
We study later the case where the variables J take the values ±N~^/'^. The two models 
coincide in the limit N ^ oo. 

We are interested in studying the statistical properties of the configurations that are 
obtained by the dynamics at zero temperature at large times. In this case the dynamics is 
such to orient the spins with the effective field hi = J2k Ji,k(^k- In other words one set 

(7i(i + l) =sign(/i,(i)). (18) 

Different algorithms differs in the order in which the rule eq. (18) is applied. 
The algorithm stops when in all the sites the following equation is satisfied 

ai{t) = sign{hi{t)). (19) 

The asymptotic configuration at large times depend on the initial configuration. We 
are interested in studied the ensemble in which the each solution of zero temperature TAP 
equation (19) is weighted with the probability of being obtained by the algorithm with a 
random choice of the initial configuration. In other words we weight each solution with 
the size of its attraction basin. 

It is known that in the large N limit the energy density of the asymptotic configuration 
does not depend on the stating configuration with probability one and it is higher that the 
ground state energy. Indeed for a sequential algorithm it is about Es = —.715 [15], while 
the ground state energy is Eq — —.7633. 
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Figure 1: Minus the average energy as function of N for the sequential, greedy and 
reluctant algoritms. The hues are linear fits. 

The simplest hypothesis would be that the set of configuration weighted with the at- 
traction basin statistically coincide with the set of all solutions of the TAP equation with 
energy equal to E. 

A precise computation of the statistical properties of the solutions of the zero tem- 
perature TAP equations for this value of the energy has not been done. However we can 
use the information we have for energies greater that EjisB = —.672 (at which an exact 
computation can be done) and for energy equal to the ground state the value. The energy 
Es is intermediate among the two so that an educated guess can be done. We could guess 
that the maximum overlap among two generic solutions ^ is is at Ejisb and it is 1 at 
Eq so that we could guess its value around .6-. 7. In general the total number of solutions 
increases about as exp{AN), where A = at Eq and A = .12 at E^sb- We guess a value 
of A in the range .04 — 0.06. These two values arc purely indicative. 

The probability distribution of the effective field P(/i) is a shifted Gaussian for energies 
greater that Erse- There are no indications that P(0) should be zero for E < Ersb- 

We will see later that some these expectations are in variance with our results coming 
from numerical simulations. We conclude that the original hypothesis is wrong and that 
the generic solution of the TAP equations, weighted with its attraction basin, is not the 
generic solution of the TAP equations of the appropriate energy. 

Before presenting the numerical results we will describe the three minimisation algo- 
rithms that we have used: the sequential algorithm, the greedy algorithm and the reluctant 
algorithm. 

• The Sequential Algorithm. 

^We will define this quantity with greater precision later on. 
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Figure 2: The average value of N{q'^) as function of N in the case of the sequential algo- 
rithm. The line is a power law fit with exponent .1; a logarithmic dependance (a straight 
hne) is also possible. 

This is the simplest algorithm to implement. One cycle of the algorithm consists in 
applying the rule (18), sequentially for increasing i from i = 1 to i = N. We repeat 
the cycles up to the moment at which a solution of the TAP equation is reached. 
This algorithm corresponds to the zero temperature hmit of an usual Monte Carlo 
or heath bath dynamics. 

• The Greedy Algorithm. 

This is the simplest algorithm to understand analytically. One step of the algorithm 
consists in applying the rule (18), for that i which minimises (7j/ij. We repeat the 
steps up to the moment at which a solution of the TAP equation is reached. This 
algorithm corresponds to the zero temperature limit of the Glauber dynamics. 

• The Reluctant Algorithm. 

This algorithm is the opposite of the greedy algorithm. One step of the algorithm 
consists in applying the rule (18), for that i which maximises aihi among those i such 
that aihi is negative. One repeats the steps up to the moment at which a solution 
of the TAP equation is reached. At each step the energy decreases, but it decreases 
of the smallest possible amount. As we shall see later this algorithm is the most 
effective in finding the configurations of smallest energies. 
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4 Numerical Results 



We have studied numerically the zero temperature dynamics of the SK model for systems 
with in the range 16-1600. We report firstly the results for sequential updating; the 
results for the other algorithms (the greedy and the reluctant) are not qualitatively different 
and they will not be discussed in details. 

For each value of many instances of the system have been generated (for 10000 
ai N — 16 to 100 at = 1600). For each choice of the coupling J we have followed 
the dynamics starting from N^^'^ different intitial configurations. For each configuration 
we have recorded the energy, the distribution of the forces hi] we have also compute the 
overlap 

q^l^^^^l (20) 

among all the iV^/^ final configurations. 

In fig. (1) we show the expectation value of the energy density E as function of N. 
The data are has been fitted as 

E = Eoo + -^, (21) 

where the values of the parameters of the fit arc £"00 = —.715, C = .25 and a = .33. The 
value of Eao differs from the gound state energy (which is -.7633). The value of a is similar 
to the one which is obtained for the N dependence of the ground state and it is compatible 
with being equal to 1/3. 

The fluctuations of the energy density go to zero when N goes to infinity approximately 
proportionally to 1/A^ (a best fit gives 1/A^'^^). 

The expectation value of q as function of goes to zero approximately as A'""'^ when 
A" — >• 00 (see fig. (2).). A behaviour of the type ^^^j^ cannot be excluded. 

These results are in good agreement with older studies [15, 16]. The surprise come from 
the study of the distribution function P{q). If we plot the function P{q) for diffeent sizes 
we do not see anything unexpected (see fig. (3)). The distribution becomes more and more 
concentrated at g = when N ^ 00. 

The crucial point is however how fast the function P{q) decreases with A^ at fixed q. 
More precisely can define the function 

/(,) = -^im!!^. (22) 

N—^00 iV 

In the ususal equilibrium SK model the function f{q) is equal to zero in the interval 
[O-^ba] and it is different from zero for q > qEA- 

The region where f{q) is zero is the essential support of the function P{q). In the 
region where /(g) =0 we can modifying the Hamiltonian by a quantity whose density goes 
to zero with A^ (e.g. imposing the appropriate bondary conditions) and we can obtain a 
system such that the value is q. 

If the function Pn{(1) becomes non trivial in the limit N ^ 00 the replica symmetry is 
broken in the usual sence. On the contrary, if the function Pn{q) becomes a delta function 
in the limit A^ going to infinity, but f{q) is zero in a finite interval, the replica symmetry 
is broken in a weak sence. 

In fig. (4,5) we plot the quantity 

rN{q)^N-^ln{['dq'PN{q')), (23) 

Jo 
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Figure 4: The results obtained with the sequential algorithm for the function r^iq) for N 
and q ranging from 16 to 1600 and from to 1 respectively. 
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Figure 5: The results obtained with the sequential algorithm for the function rjv(?) for N 
and q ranging from 100 to 1600 and from to .4 respectively. 



with u = 2/3. 

This quantity has a very weak dependence on N and it is rather likely that it goes 
to a non zero value when N goes to infinity in the whole interval 0-1. If we accept this 
conclusion, we find that PNi^) decreases slower than an exponential for all values of q, 
included q — 1 and therefore the function /(g) is always zero in the whole interval 0-1. 
Replica symmetry is thus broken in a weak sence. 

If rjv(l) has a finite limit, as it is suggested by this data, the probability of finding 
two times the same solution goes to zero slower than an exponential of A^, implying that 
the effective numer of different solutions which can be reached with this method diverges 
slower that an exponential of A^. 

The other quantity for which we obtain surprising results is the probabily distribution 
of the fields hi = J2kJi,kO'k- The function Piv(^) depends weakly on N (the expectation 
value of h is the energy. 

The function P]s[{h) seem to vanish linearly at /i = 0. In order to show this behaviour 
in fig. (6) we plot the function 

for = 1600 where the constant C (which vanishes when — ^ oo) has been choosen in 
such a way to have a smooth behaviour at small h. We have found that C = 2A^^^/^ is a 
good choice. This value of C is natural, indeed the variable h make take values (n-|- 1/2)C, 
with n integer. 

A linear behaviour of the function P{h) is similar to the one observed at thermal 
equilibrium at zero temperature, where P{h) is proportional to h. 
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Figure 6: The function P{h)) in the case of the sequential algorithm. 



We do not have a simple explanation for such behaviour. The process of minimizing 
the energy of a given spin has the side effect of decreasing also the energy of the other 
spins, pushing the distribution of h far from zero. Howevever it is unclear how to transfom 
this observation in a quantitive prediction. 

We have anahzed also the other dynamics ad we have obtained comparable results. The 
only remarkable effect that the dependence of the energy on the algorithm is counterintu- 
itive. Indeed the greedy algorithm is the worst (the energy is higher of about .007 than 
the sequential algorithm) and the reluctant algorithm is much better (the energy -.746 is 
lower of about .03 than the sequential algorithm). This effect is likely due to the fact that 
the greedy algorithm is more easily trapped in local minima with high energies and that 
the slower reluctant algorithm avoids these traps. 

5 Conclusions 

The probability distribution of the configurations obtained by the zero temperature dy- 
namics have rather peculiar statistical properties. The results obtained are very different 
from those coming taking a generic solution of the TAP equations with the appropriate 
energy. The maximum allowed value of q (from the thermodynamic point) is 1, not smaller 
than 1 as it would happen for the generic solutions of that energy. Morover the distribution 
of h is zero at h =0, while there are no reasons that P{h) should be equal to zero at h — 
in the case if the generic solution. 

It is rather hkely the corresponding effective Hamiltonian is not so simple and its 
statistical properties should be computed using a replica symmetry breaking approach. 

It would also interesting to connect the small h behaviour with the power laws observed 
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in the dependence of physical quantities (hke the energy of the remanent magnetization) 
or the number of iterations in the sequential dynamics. 
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